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THE SYMBOLIC MANIPULATION LANGUAGE FORMAC 
by Paul Swlgert 
Lewis Research Center 

SUMMARY 

A method for constructing quadrature coefficients through the use of a FORMAC com- 
puter program is presented. The computer program has the ability to construct exact 
rational fraction coefficients for quadrature formulas of a wide class of weight functions. 
Quadrature coefficients, constructed by the FORMAC program, are presented for certain 
weight functions and collocation points. 


INTRODUCTION 

Numerical integration of functions is usually easy and straightforward on high-speed 
computers. This ease, obtained through the use of a general purpose integration sub- 
routine, may be misleading as to efficient programing. The general subroutine must ap- 
ply an unweighted quadrature formula, such as Simpson’s rule, many times over the inte- 
gration interval to achieve a solution. To keep the error in the solution small some 
method of testing the accuracy, like differencing, and halving the increment must be in- 
corporated. The very generality of such a program leads to an increase in computer 
time needed to perform a specific integration. This is especially true if the integrands 
are not well behaved (e. g. , l/y/x and lnx). Functions such as these cannot be approxi- 
mated near the origin by a polynomial and thus cause excessive computer time when being 
integrated by an unweighted quadrature formula. 

A more efficient method of performing numerical integration may be done by quadra- 
ture formulas developed specifically for the function to be integrated. The efficiency of 
a specific quadrature formula over a general purpose integration subroutine increases as 
the number of integrations to be performed goes up. Thus much computer time may be 
saved if a single quadrature formula can be used for integrands of a common functional 
form (e. g. , f (x)yjx). 



H=(Hi H 2 Hg • • • H n ) 

It can be shown that, if all x i are independent, the inverse of the Vandermonde matrix 
may be computed from the following equations: 

V = LU 


where U and L are upper and lower triangular matrices, respectively. The elements 
Uy of U~* are 


u ij 

i > j 

“ij = 1 

i = j 

u l,] = - u l,i-l x i-l 

1 < j 


i < i 


The elements f y of L~* are 

*13 = 0 1 < i < j 


£ 11 “ 1 


k=l 

k*j 


i > j > 1 


k=l x j " x k 


The final form of the desired matrix equation may now be written 

HF T = MU" 1 L“ 1 F T 


The right side of this equation is computed by the FORMAC program described in the 
following section. 


4 



DESCRIPTION OF FORMAC PROGRAM 


The function of the FORMAC program is to construct quadrature formulas based on 
the method described in the preceding section. This computer program, consisting of a 
main program and an integration subroutine, is written in FORMAC (Formula Manipula- 
tion Compiler) for the IBM 7094-7044 direct couple system. FORMAC was chosen as the 
programing language because of its ability to operate in rational fraction arithmetic and 
its capability to manipulate mathematical expressions. The former ability allows the 
calculation of exact quadrature coefficients, while the latter allows for simple input and 
closed form integration of the various weight functions. 

Listings and flow charts of the main program (INTCOF) and the integration sub- 
routine (POLINT) may be found in the appendix. The general flow of the main program 
and integration subroutine may be obtained from the flow charts. All the matrices and 
vectors necessary are constructed and manipulated in the main program, while closed 
form integration, used in construction of the moments of the weight function, is per- 
formed in the integration subroutine. 

The integration subroutine will integrate, in closed form, any linear combination of 
the two following functions: 


r 

gw = 


k=l 



h(x) = ^2 C k x k In x (d k * -1) 
k=l 


where C k and d k are rational constants and r and s are positive integers. This 
subroutine relies on the following indefinite integral formulas: 



(p/q)+l 
dx =- 

P+ 1 


q 


(p/q* -i) 



dx = In x 
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In x dx 


-(p/q)+l 

k 

P+1 

q 


In x - 


x (p/q)+l 



(p/q * -i) 


The result of integration by this subroutine is a symbolic expression for the integral. 
This integral must, however, be evaluated at the upper and lower limits to determine the 
moments. This evaluation is performed in the main program. 

Since lim x In x = 0, it can be seen that, if the integration range is (0, 1) and 
x— 0 

p/q > -1, irrational numbers will not have to be calculated in any of the integrations. 

For any other range of integration, it is necessary to keep such quantities as 2^^ and 
In 3, obtained from applying the limits of integration in determining the moments, in 
their symbolic notation so that exact coefficients may be obtained. This is done auto- 
matically in the main program. 


PROGRAM USE 

Input to the quadrature program is entered on two cards. The first card contains 
the number of collocation points n. This number is in fixed point notation and must end 
in card column 10. The other input goes in the first 72 card columns of the second card, 
the order and form of which is shown in table I. Note that the weight function is in terms 

TABLE I. - DEFINITION OF INPUT 


Notation 

Description 

Example 1 

Example 2 

w(x) 

Weight function 

Z**(-2/3)$ 

Z**3*FMCLOG(Z)$ 

a 

Lower limit of 
integration 

0$ 

0$ 

b 

Upper limit of 
integration 

1$ 

1$ 

x i 

Collocation 

points 

l/4$l/2$3/4$l$ 

0$l/4$l/2$3/4$l $ 


of a dummy variable Z, that FMCLOG(Z) is the FORMAC notation for In z, and that 
each "piece” of this input is terminated by a dollar sign. Blanks are ignored and may be 
used freely on this card. The restriction to the first 72 card columns may be altered by 
changing the READ and DIMENSION statements in the main program. The two examples 
from table I as they would appear on cards, along with their corresponding output from 
the program are shown hereinafter. 
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4 

Z**<-2/3)$ OS 1 $ 1 / 4 $ i/2S 3/4$ 1$ 

5 

Z**3*FMCL0G(Z ) $ o$ IS OS 1/4$ 1/2$ 3/4$ 1$ 


THE WEIGHT FUNCTION IS Z**tI-~2)*3**(-l))$ 
THE LOWER INTEGRATION LIMIT IS 0$ 

THE UPPER INTEGRATION LIMIT IS IS 
THE 'X* POINTS FOLLOW 


4**(-l) $ 

2** (- l ) S 
3*4** (-1 )$ 

IS 

THE 4 POINT INTEGRATION 0 UAUKA TORE FUR THE ABOVE CONDITIONS IS - 


213*3 5**<-l)*YI 1 )- 5 1 3* 70** (-1) *Y 12 ) +27*5** (-1 ) *Y ( 3 )-8 1 * 70 ** ( - I } * Y ( 4 ) S 


THE WEIGHT FUNCTION IS Z**3*FMCL0G ( Z ) S 
THE LOWER INTEGRATION LIMIT IS OS 
THE UPPER INTEGRATION LIMIT IS IS 
THE 'X' POINTS FOLLOW 


0 $ 

4* * ( — 1 ) $ 

2** (-1 ) $ 

3*4** (- 1 ) S 
IS 

THE 5 POINT INTEGRATION QUADRATURE FOR THE ABOVE CONDITIONS IS - 


Y 1 I)*5*2Il6ti**t-II+YI2I*|-21B) *3 3 0 75 ** (-1 ) + Y I 3 > * ( -208 I * 1 1025** 1 - I I +Y 141*1-1581*4 72 5** l-l) + YI5)*( 
— 12 ( I*33075**t — ITS 


Certain restrictions are inherent in this program. The specified integration interval 
and collocation points must be entered as rational fractions. Only weight functions of the 
type x P//q or x p//q In x (p/q > -1) are acceptable to the program. Since this program 
operates in rational arithmetic, the possibility of integer overflow is great. To keep this 
possibility at a minimum, the integration range should be (0, 1). No more than nine col- 
location points should be used, and collocation points with large numerators and denomi- 
nators should be avoided. 
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DISCUSSION 


Quadrature coefficients that have been constructed by this method are presented in 
table n (pp. 19 and 20). Table n(a) lists the coefficients, fL, for equally spaced colloca- 
tion points over the closed integration interval [0, 1] and the weight functions 1, x'*'/ 2 , 
x -1 / 2 , In x, and x In x. It might be noted that a linear transformation is all that is 
needed to transform any finite interval, on the real line, to the interval [0, 1]. The coef- 
ficients presented, therefore, are applicable to any finite integration interval. The 
FORMAC program may be used, however, to generate quadrature formulas for other in- 
tervals and weight functions. 

To give an example of the use of such a linear transformation and to show the use of 
a quadrature formula from table 11(a), suppose a function F(4) is to be calculated and is 
given by 


F(4) = 



_ — * 

8 - o 1/2 


where 4 is finite and f(t) is not singular in [0,4]. From the transformation 


t = 4 (1 - z) 


the following equations are obtained: 


dt = -4 dz 


1 


z =< 


0 


when t = 0 
when t = 4 


Substituting for t and dt 


gives 


F(4) 


-5 



f[4 (1 - z )]dz 
(4 - 4 +£z) 1/2 


8 


and, finally, 


F(!) = « 


-*l/2 


/ 

y 0 


z -1 / 2 f[£(l - z)]dz 


This integral may be readily evaluated by using a quadrature formula with a weight func- 
tion of x -1 / 2 from table H(a). 

Quadrature coefficients for the same weight functions as those in table n(a) are given 
in table n(b) for equally spaced points over the half open integration interval (0, 1]. 

This table is useful for integrating functions that have a numerical singularity at the ori- 
gin. For example, consider the numerical integration of the following: 


/ 


V si 


sin x 


dx 


This may be rewritten as 


/ I /sin x 
yfc 1 x 


dx 


At x = 0, (sin x)/x cannot be evaluated on a computer. To handle such an integration, 
a quadrature formula with a weight function of x -1 / 2 may be chosen from table n(b). 

By this method the need to evaluate (sin x)/x at x = 0 is eliminated. 

The two examples given previously are only presented to show how the tables can be 
used. There are, of course, many more applications of the quadrature formulas pre- 
sented. 

In integrating some functions, it may be desirable to place more collocation points 
in a region of rapid functional change than in other regions. This may be accomplished 
by specifying the desired points in the program and obtaining a quadrature formula that 
is specifically designed for the function. This method might also be applied to a function 
given in tabular form but at nonequidistant points. 

Errors resulting from truncation and round off associated with floating point arith- 
metic have been eliminated from the calculation of these coefficients by the use of 
FORMAC. To check the errors that can accumulate from floating point round off and 
truncation, the presented program was run using single precision floating point arith- 
metic. hi five-point quadrature formulas, the errors in the coefficients ranged from 
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one in the fifth to one in the eighth significant figure. The accuracy could be improved 
by making use of double precision and better matrix inversion techniques; but the diffi- 
cult task of determining exactly what the error is in the coefficients would still remain. 

A discussion of errors occurring in approximating the integral of a function by a 
quadrature formula may be found in reference 4. Reference 2 gives exact quadrature co- 
efficients for the weight function w(x) = f or n _ 2(1)11 along with an excellent 

discussion of error bounds on quadrature formulas. 


CONCLUDING REMARKS 

The calculation of exact quadrature coefficients has been reduced to a FORMAC com- 
puter program. This program may easily be used to produce quadrature formulas for a 
large class of weight functions. The use of quadrature formulas designed for a particular 
integrand promises to save computer time over the use of a general purpose integration 
program. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 29, 1966. 
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APPENDIX - FORMAC PROGRAM LISTINGS 


SIBFMC INTCRF 

FORMAC PROGRAM TO EVALUATE N POINT INTEGRATION 
COEFFICIENTS USING THE VANDERMONDE MATRIX 
SYMARG 

DIMENSION IX(12»12)»IY(12)»JY(12)»KY(12)»X(12)»IN(12) .OUT( 17) 
ATOMIC Y< 12) ,Z,LN1,LN2 
LOGICAL KODE 

INTEGER QUAD»X,Y»W,Z»X1 ,X2»T1 ,T2 
FMCDMP LATER 
C READ INPUT 

1 READ (5.101) N. ( INI I ) ,1 = 1 .12 ) 

101 FORMAT (110/(1 2 A6 > > 

J = 0 

LET W = ALGCON IN ( 1 ) » J 
LET XI = ALGCON I N ( 1 ) ,J 
LET X2 = ALGCON INI 1 ) ,J 
DO 2 1=1, N 
I = I 

LET X ( I ) = ALGCON I N ( 1 ) ,J 

2 CONTINUE 

C CONSTRUCT MOMENT VECTOR 

DO 10 1=1, N 
I = I 

LET T 1 = W*Z**(I-1) 

CALL POLINT ( T1 »Z ,KODE) 

IF ( .NOT. KODE ) STOP 

LET KY(I> = SUBST Tl.(FMCLOGIZ) ,LN2) »(Z,X2) 

LET T 1 = SUBST T 1 , ( FMCLOG ( Z ) , LN 1 ) , ( Z . X 1 ) 

LET KY ( I ) = FXPAND KY(I> - T1 

10 continue 
ERASE T 1 

C CONSTRUCT LOWER INVERSE 

DO 30 I = 1 ♦ N 
I = I 

DO 30 J= 1 , N 
J = J 

LET I X ( I , J ) = 1 
I F ( I.GE.J) GO TO 19 
LET I X ( I , J ) = 0 
GO TO 30 

19 DO 20 K = 1 , 1 
K = K 

IF(K.EQ.J) GO TO 20 

LET I X ( I , J ) = EXPAND I X ( I , J ) * ( 1 / ( X ( J ) -X ( K ) ) ) 

20 CONTINUE 
30 CONTINUE 

C MULTIPLY COLUMN VECTOR BY LOWER INVERSE 

DO 40 1 = 1, N 
I = I 

LET I Y ( I ) = 0 
DO 40 J=1 »N 
J = J 

LET I Y ( I ) = EXPAND I Y ( I ) + IX(I»J)*Y(J) 

ERASE I X ( I . J ) 

40 CONTINUE 

C CONSTRUCT UPPER INVERSE 

LET I X ( 1 , 1 ) = 1 
DO 50 J=2,N 
J = J 

LET IX(1,J) = EXPAND IX( 1 ,J-1 >*X( J-l )*(-l > 

50 CONTINUE 
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I X ( I , J-l ) * X ( J-l ) 


DO 60 1=2 ,N 
I = I 

DO 60 J=1 »N 
J = J 

I F ( I.LT.J) GO TO 59 
I F t I .EQ.J) GO TO 58 
LET I X ( I * J ) = 0 

GO TO 60 

58 LET IX l I » J ) = 1 
GO TO 60 

59 LET IX(I»J) = EXPAND IX(I-l.J-l) - 

60 CONTINUE 

C MULTIPLY COLUMN VECTOR BY UPPER INVERSE 

DO 70 I =1 »N 
I = I 

LET JY ( I ) = 0 

DO 70 J=1 *N 
J = J 

LET JY ( I ) = EXPAND JY(I> + IX(I.J) *IY(J) 

ERASE IX(I,J) 

7 C CONTINUE 

DO 80 1=1, N 
I = I 

ERASE I Y ( I ) 

80 CONTINUE 

C MULTIPLY RESULTING COLUMN VECTOR BY THE MOMENT VECTOR 

LET QUAD = 0 
DO 90 1 = 1, N 
I = I 

LET QUAD = QUAD + KY(I)*JY<I> 

ERASE KY ( I ) , JY ( I ) 

90 CONTINUE 

LET QUAD = EXPAND QUAD 

LET T 1 = QUAD 

ERASE QUAD 

LET QUAD = 0 

DO 91 1 = 1, N 

I = I 

LET T 2 = COEFF T1 * Y ( I ) 

LET QUAD = QUAD + Y(I)*T2 

91 CONTINUE 
FRASE T1 » T? 

LET QUAD = ORDER QUAD, INC, FUL 

LET KODE = FIND QUAD ,APP .ONE ,( LN 1 , LN2 ) 

I F ( .NOT .KODE ) GO TO 93 
LET KODE = MATCH ID.X1.1 
I F ( .NOT. KODE) GO TO 92 
LET QUAD = SUBST QUAD,(LN1,0) 

AUT SIM ON 

LET QUAD = SUBST QUAD , ( LN2 , FMCLOG ( X2 ) ) 

GO TO 93 

92 LET KODE = MATCH ID.X2.1 I 

IF( .NOT. KODE) GO TO 94 

LET QUAD = SUBST QUAD,(LN2,0) 

LET KODE = MATCH ID.X1.0 
IF(KCDE) GO TO 93 
AUTSIM ON 

LET QUAD = SUBST QUAD ,( LN1 , FMCLOG ( X 1 ) ) 

GO TO 93 
94 AUTSIM ON 

LET QUAD = SUBST QUAD ,( LN 1 , FMCLOG ( X 1 )).( LN2 .FMCLOG ( X2 ) ) 
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C WRITE OUTPUT 

93 QUEST = 0.0 

LET QUEST = BCDCON W, OUT, 17 
WRITE (6,201) (OUT< I) ,1=2*17) 

QUEST = 0.0 

LET QUEST = BCDCON XI, OUT, 17 
WRITE (6,202) ( OUT ( I ) , 1 = 2 ,17 ) 

QUEST = 0.0 

LET QUEST = BCDCON X2,0UT,17 
WRITE (6,203) ( OUT ( I ) , I =2 , 1 7 ) 

WRITE (6,204) 

DO 99 1=1, N 
I = I 

QUEST = 0.0 

LET QUEST = BCDCON X( I), OUT, 17 
WRITE (6,205) ( OUT ( J ) , J=2 , 17 ) 

99 CONTINUE 

WRITE (6,206) N 
QUEST = 0.0 

100 LET QUEST = BCDCON QUAD, OUT, 17 
WRITE (6,207) ( OUT ( I ) , I =2 ♦ 17 ) 

IF(QUEST.NE.O.O) GO TO 100 
ERASE QUAD 
AUTSIM QNINT 
GO TO 1 

201 FORMAT ( 24H1 THE WEIGHT FUNCTION IS ,16A6) 

202 FORMAT ( 32HKTHE LOWER INTEGRATION LIMIT IS »16A6) 

203 FORMAT ( 32HKTHE UPPER INTEGRATION LIMIT IS »16A6) 

204 FORMAT ( 22HKTHE -X- POINTS FOLLOW,///) 

205 FORMAT (25X,16A6> 

206 FDRMAT ( 5HKTHE ,I2,59H POINT INTEGRATION QUADRATURE FOR THE ABOVE 
♦CONDITIONS IS -,///) 

207 FORMAT < 1HK , 10X , 16A6 ) 

END 
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Convert BCD 
input to FORMAC 
internal repre- 
sentation 


Construct 

moment 


Construct 

lower 

inverse 


Multiply 
column vector 
by lower inverse 


Construct 

upper 


f fU 7U ! 

Replace all 

Multiply column 

Multiply resulting 

FMCLOG(l) 

vector by upper 

► column vector by — ► 

by 0 in re- 

inverse 

moment vector 

sulting quad- 



rature 

Figure 1. 

- Flow chart for main program. 










SIBFMC POLINT 

SUBROUTINE POL I NT C Y ,X , KODE ) 

C FORMAC SUBROUTINE FOR CERTAIN CLOSED FORM INTEGRATIONS 

INTEGER X*Y » Y 1 * Y2 *Y3 » Y4 » D I V 
LOGICAL QUEST, KODE 
SYMARG Y,X 
ATOMIC X 

LET Y = EXPAND X*Y 
LET N = CENSUS Y » TERM 
LET Y 1 = 

DO 10 1*1, N 
K = 0 

I F ( I .EO.N) GO TO 1 
LET Y 2 = PART Y,M 
IF(M.EQ.4> GO TO 2 
1001 FMCDMP NOW 

KODE = .FALSE. 

RETURN 

1 LET Y2 = Y 

2 LET QUEST = FIND Y2 , APP , ALL * ( FMCLOG ( X ) ) 

IF ( .NOT. QUEST) GO TO 3 

LET Y 2 = EXPAND Y2/FMCL0GIX) 

K = 1 

3 LET Y3 * Y2 

4 LET Y4 = PART Y3*M 
IF(M.NE.l) GO TO 5 
LET Y 3 = Y4 

GO TO 4 

5 IFtM.NE.5) GO TO 6 

LET QUEST = FIND Y3 ,APP , ALL . ( X ) 

IF(QUFST) GO TO 4 
LET Y3 = Y4 

LET QUEST = FIND Y3 ,APP , ALL * ( X ) 

IF(QUEST) GO TO 4 
IF(K.EO.l) GO TO 1001 

LET DIV = Y2*FMCL0G ( X ) 

GO TO 9 

6 IFIM.NE.2) GO TO 7 

IF(K.EQ.O) GO TO 61 

LET DIV = Y2*(FMCL0G(X)-1 ) 

GO TO 9 

61 LET DIV = Y2 
GO TO 9 

7 IFfM.NE.3) GO TO 8 

IF(K.EQ.l) GO TO 1001 

LET DIV = Y2*FMCL0G t X ) 

GO TO 9 

8 IF(K.EQ.O) GO TO 81 

LET DIV = Y2*(FMCL0G(X)/Y3— 1/Y3**2) 

GO TO 9 

81 LET DIV = Y2/Y3 

9 LET Y 1 * Y 1 + DIV 
10 CONTINUE 

LET Y = EXPAND Yl 
ERASE Y1*Y2*Y3*Y4»D1V 
KODE = .TRUE. 

RETURN 

END 


From cah vector 
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TABLE H. - QUADRATURE COEFFICIENTS Hj 
(a) Closed interval, [0, 1]; x. = (i - l)/(n - 1) (i = 1, 2, ■ • • , n) 


Weight Number of Functional values Weight Number of Functional values 


function, 

w(x) 

collocation 

points, 

n 

f(Xj) 

I(x 2 ) 

J(x 3 ) 

f(x 4 ) 

f < x 5> 

f(x 6 ) 

f(x ? ) 

«X 8 ) 

« x 9> 

function, 

w(x) 

collocation 

points, 

n 

f(x x ) 

f( x 2 ) 

f< x 3 ) 

f(x 4 ) 

f ( x 5? 

f(x 6 ) 

f(x 7 ) 

f (x 8 ) 

« x 9 ) 

1 

2 

1 

1 








In x 

2 

^3 

^1 










2 

2 










4 

4 









3 

1 

2 

1 








3 

-17 

-5 

_ 1 _ 









6 

3 

6 









36 

9 

36 








4 

1 

3 

3 

1 







4 

-11 

-19 

-1 

-1 








8 

8 

8 

8 








32 

32 

32 

32 







5 

1 _ 

16 

_2 

16 

_7_ 






5 

-727 

-398 

-1 

-98 

23 







90 

45 

15 

45 

90 







2700 

675 

225 

675 

2700 






6 

19 

25 

25 

25 

25 

19 





6 

-767 

-1975 

25 

-425 

25 

-13 






288 

96 

144 

144 

96 

288 






3456 

3456 

576 

1728 

3456 

1152 





7 

41 

JL 

9 

34 

9 

_9 

41 




7 

-2659 

-2721 

4293 

-8083 

2823 

-369 

1439 





840 

35 

280 

105 

280 

35 

840 





14112 

4900 

39200 

22050 

39200 

4900 

352800 




8 

751 

3577 

49 

2989 

2989 

49 

3577 

751 















17280 

17280 

640 

17280 

17280 

640 

17280 

17280 


x In x 

2 

-5 

J. 









9 

989 

2944 

-464 

5248 

-454 

5248 

-464 

2944 

989 



36 

9 










28350 

14175 

14175 

14175 

2835 

14175 

14175 

14175 

28350 


3 

^1 

24 

-7 

36 

-1 

72 







X i/2 

2 

_4_ 

2 









4 

-77 

-107 

-17 

-71 








15 

5 










3600 

800 

200 

7200 







3 

4 

16 

_6_ 








5 

-13 

-67 

-67 

-41 

-19 







105 

35 

35 









1080 

675 

900 

675 

5400 






4 

4 

_6 

12 








6 

-1385 

-6215 

-5455 

-2785 

-5945 

-263 






105 

35 

35 

35 








169344 

84672 

84672 

42336 

169344 

84672 





5 

4 

64 

368 

448 

274 






7 

-47 

-289 

-1859 

-787 

-44 

-143 

-107 





297 

385 

3465 

1485 

3465 







8400 

4900 

39200 

11025 

1225 

4900 

70560 




6 

2036 

335 

4280 

2480 

620 

8429 

















135135 

3861 

27027 

27027 

2457 

135135 
















7 

76 

2384 

178 

96 

188 

5776 

1236 
















10725 

25025 

5005 

455 

5005 

25025 

25025 















v-1/2 


4 

3 

4 

5 

68 

105 

100 

189 

36 

77 

1476 

3575 


2 

3 

16 

15 

6 

7 

832 

945 

1685 


2_ 

15 
12 
35 

16 
315 

40 


16 

105 

64 

135 

40 


62 

945 

460 


23 


20688 

-1062 

10592 

-516 

8208 

3176 




25025 

5005 

15015 

5005 

25025 

75075 
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